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ABSTRACT 


In the source-free mant le/ frozen- f lux core magnetic earth model, 
the non-linear inverse steady motional induction problem has been solved 
using the method presented in Part IB. The present paper describes how 
that method has been applied to estimate steady, broad-scale fluid 
velocity fields near the top of Earth's core that induce the secular 
change indicated by the Definitive Geomagnetic Reference Field (DGRF) 
models from 1945 to 1980. Special attention is given to the derivation 
of weight matrices for the DGRF models because the weights determine the 
apparent significance of the residual secular change. The derived 
weight matrices also enable estimation of the secular change signal-to- 
noise ratio characterizing the DGRF models. 

Two types of weights were derived in 1987-88: radial field weights 
for fitting the evolution of the broad-scale portion of the radial 
geomagnetic field component at Earth's surface implied by the DGRFs , and 
general weights for fitting the evolution of the broad-scale portion of 
the scalar potential specified by these models. The difference is non- 
trivial because not all the geomagnetic data represented by the DGRFs 
constrain the radial field component. For radial field weights (or 
general weights), a quantitatively acceptable explication of broad-scale 
secular change relative to the 1980 Magsat epoch must account for 
99.94271% (or 99.98784%) of the total weighted variance accumulated 
therein. Tolerable normalized root -mean-square weighted residuals of 
2.394% (or 1.103%) are less than the 7% errors expected in the source- 
free mant le/ frozen- f lux core approximation. 
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1. INTRODUCTION 


The work of ELSASSER (1946a, b 1947 ) and others has led to general 
acceptance of the importance of induction effects in geomagnetism (see, 
e.g., JACOBS, 1987). The importance of steady induction effects in 
geomagnetism, although clearly less, cannot be dismissed simply because 
slow changes of the geomagnetic field, the secular variation (SV) , are 
not purely linear with time— as evidenced by the inaccuracy of 
predictions by BOND (1668) . Unlike steady magnetic flux diffusion, 
motional induction by steady fluid flow near the top of the electrically 
conducting liquid outer core need not imply constant SV. Consider, for 
example, a high conductivity core threaded by a partly non-axisymmetric 
magnetic field of internal origin: steady, uniform westward flow of this 
core relative to a magnetically source-free mantle induces unsteady SV 
in the mantle reference frame. 

1 . 1 Synopsis 

Part IA (VOORHIES, 1992) offered some reasons to consider steady 
induction effects in geomagnetism. The source-free mant le/ frozen- f lux 
core (SFM/FFC) earth model (wherein a rigid, impenetrable, electrically 
insulating mantle of uniform magnetic permeability surrounds a 
spherical, inviscid, and perfectly conducting core in anelastic flow) 
was introduced and used to focus attention upon the fluid motion near 
the top of the core. The theory underlying some estimates of core 
surface flow was summarized. Some implications of a few kinematic and 
dynamic hypotheses were derived. The solution to the forward steady 
motional induction problem was found to be an example of deterministic 
chaos. Special attention was given to the long-term geomagnetic effects 
of persistent, surficially geostrophic flow. 

To investigate steady induction effects in a quantitative way, it 
is useful to treat the supposition of steady surficial core flow as if 
it were a hypothesis to be tested geomagnet ical ly . Then the idea is to 
see if such a flow can be designed to induce the secular change signal 
indicated by geomagnetic data (or broad-scale spherical harmonic models 
thereof) and, if not, assess the significance of the residual, unmodeled 
signal. Part IB (VOORHIES, 1993) offered a fairly general method for 
estimating a steady fluid flow at the top of the core in the context of 
the SFM/FFC model. The method is based upon iterative solution of the 
linearized weighted least -squares problem, but admits optional biases 
favoring surficially geostrophic flow and/or spatially simple flow. 

As noted in Part IB, solution of the non-linear inverse problem 
posed by the hypothesis of (piecewise, statistically) steady core 
surface flow requires specification of a complete initial geomagnetic 
condition: the radial field component everywhere at the top of a FFC, 
B r (b,t), or, for a SFM, at Earth's surface, B r (a,t). Such a condition 
will not be specified by discrete geomagnetic data alone within the 
forseeable future. Moreover, tests of steady flow are likely embedded 
in some magnetic earth model-be it the SFM/FFC model or a refinement 
thereof. The rigor of a purely geomagnetic test is thus limited: it may 
be that errors in the initial condition or in the earth model, rather 
than errors in the steady motions hypothesis, lead to unacceptably large 
residuals. It is thus doubtful that geomagnetic tests of the steady 
motions hypothesis can be absolutely decisive. Fortunately, we already 
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know that both the SFM/FFC model and the steady motions hypothesis are 
but conceivably relevant idealizations. The question remains as to what 
fraction of recent SV cannot be readily attributed to steady induction 
effects. Such unmodeled signal may contain evidence for many other 
interesting geomagnetic effects. It turns out that application of the 
method developed in IB is well suited to answering this question. 

1.2 Review of the Theory 

For steady flow, the radial component of the magnetic induction 
equation at the top of a frozen-flux core of radius b is, in spherical 
polar coordinates (r,0,<|>), 

a t B rp (b,t) + v s (b)*V s B rp (b, t) = B rp (b, t)3 r u(b) (0) 

which is a special case of the ROBERTS & SCOTT (1965) equation. With 
sums running over repeated subscripts, the compact spherical harmonic 
expansions for the input radial magnetic flux density B r (r,t) and for 
the predicted radial magnetic flux density B rp (r,t) are, respectively, 

B r (r,t) = g i (r,t)S i (0,<|>) < la > 

B rp (r,t) = Y i (r,t)S i (0,(|)) . (lb) 

Here S i (0,<|)) is a Schmidt-normalized spherical harmonic; compact index i 
is uniquely determined by, and uniquely determines, both the spherical 
harmonic degree n(i) and the order m(i) . Radii a = 6.3712 Mm and b = 
3.480 Mm are of interest here, so define: g^ = g^(a,t); = g^(b,t); y ^ 

= Yi (a, t) ; and T ^ = Yi (b, t) . For a SFM, g^_ = ^ij G j and Yf = * 

where Y — = [ b/a ] n ^ + 2 8 i j is an element of the diagonal upward 
continuation matrix. The g i are the input Schmidt -norma 1 i zed Gauss 
coefficients (g n m ,h n m ) for the usual internal scalar magnetic potential 
at reference radius a multiplied by [ n ( i ) + 1]. 

The hypothetically steady fluid flow at the top of the core 

A 

v s (b) = V s T(b) x r + V s U(b) (2a) 

is parameterized in terms of the spherical harmonic coefficients for the 
expansions of the streamf unction -T and velocity potential -U: 

T(b) = « i S i (0, <(>) U(b) = (0,^) - (2b) 

Equation (IB. 15) follows from substitution of (la-2b) into (0): 

^t^k = ^i X ijk^ a j + ^i Y ijk^Pj = P kj a j + ^kjPj • 

The (X i j k ,Y i j k ) are given by (IB. 16); (PfcjiQfcj) are given by (IB. 17). 
Upward continuation of (3) gives equations (IB.l8a-c): 


3 t Yi - “ hk^ ^i x ijk* a j + *^i Y ijk^Pj^ 

(4a) 

= hkt p kj a j + ®kjPj^ 

(4b) 

= Pijttj + qijPj = A n ^. 

(4c) 
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The matrix notation for (4c) is d t y = A^, where is the extended 

column vector of streamfunction and velocity potential coefficients and 
A is the matrix of normal equations coefficients. The elements A-, of A 
vary with time due to their dependence upon B (b,t), hence upon % and 
the initial geomagnetic condition. ~ 

1.3 Summary of the Method 

The approach to estimating hence the a^, the (}-, T, U, and 

steady flow v g (b), described in Part IB is based on minimizing the 
function 

4*A 2 = 4K(A r 2 + A.gAg 2 + X^) . (5a) 

In (5a), Xg and are adjustable damping parameters. The semi- 
normalized square-weighted residual, or misfit, is A r 2 : 

t £ ip 

4jcA r 2 = I (g - y) W ( g - y)dt 


t f t t 

= 1 { [J 0 T g - A^)dT] T W[J (d x g - A^)dT]}dt 

fc o fc o = fc o = ' 


35 [3 T g - A^] T W[3 T g - A^] 


(5b) 


where W is the time-dependent radial field weight matrix and the initial 
condition is taken to be g(a,t Q ) = Y(a,t 0 ). For general weights, £2// 
replaces W/4 n in (5b). The mean-square ageostrophy of the flow is A g 2 T 

2 n n v(b) 

4jtX g Ag 2 = A.g J J (d r u(b)cos0 + s i n8 ) 2 s i n0d0d<|> 

0 0 b 

v 

= 4rtXg<{3 r ucos0 + -sin0) 2 > 

b 

= [B4) T A g [B^]. (5 C ) 

The mean-square radial vorticity plus the mean-square downwelling is 
just A d z = < [« r (b) ] 2 + [3 r u(b)] 2 >: 

47tA. d A d 2 = 4 T A d ^. ( 5d) 


The linearized weighted least -squares estimate of the parameters is 
obtained by treating A as if independent of | and minimizing (5a) : 
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( 6 ) 


£ = [A t W A + B T A g B + A d J 1 [A t W 3 x g] 


which is (IB. 32) . 

Because A depends on the predicted radial field component B r p(b,t), 
it depends upon I\ hence upon \ and the initial condition. Minimization 
of (5a) is thus a non-linear inverse problem. To solve this problem by 
iterative solution of the linearized problem, initial matrices A(0) are 
computed from the input g(t). Then (6) is solved for £(1) which, in 
turn, is used to construct v g (b) and solve the forward problem (0) from 
the inital condition at t Q . The resulting predictions B r p(b,t) are used 
to compute the time-varying predicted radial field coefficients 7(1) . 
the residuals 6g(l) = g - Y(l), and new matrices Ad) for the next 
estimate, \ ( 2 ) . ~ With 6g_ ( j ) = g-Yjj> and £(j + l) = | ( j ) + 5^( j + 1 ) , 
the correction to |(j) is given by deep iteration equation (IB. 33a): 


5£( j + 1) = [A(j) T W A ( j ) + B T AgB + A d ] ^ [ A ( j ) *^ W 8d x g(j)] 
- [ B^ Ag B] [% ( j ) - ^g Q ] _ [^ ( j ) “ ^do^ ' 


(7) 


Shallow iteration is described by equation (IB. 33b) . 

Note that % o is taken to be 0 in the calculations. The routine 
prior bias Is also taken to be 0; however, some calculations have 

been performed with a learning algorithm so as to escape from this bias 
against non-trivial parameters. Then | do is taken to be |(j) and 
A., — * A. d (j + l) is selected to keep ^(j + 1) close to ^(j) (see IB, section 
3 . 3 ) The learning algorithm is used to explore the tightness of fit 
afforded by steady flows which appear so spatially complicated as to 
inhibit convergence of the iteration scheme (7)— as expected if (5a) has 
multiple extrema for very small values of X d - 

2. APPLICATION 


The foregoing method has been applied to estimate piecewise steady, 
optionally geostrophic, and optionally damped, surficial core motions 
which induce the secular change indicated by the Definitive Geomagnetic 
Reference Field (DGRF) models. The DGRF models were derived by Working 
Group 1.1 of the International Association of Geomagnetism and Aeronomy 
( I AG A, 1988). The DGRFs are slowly varying, broad-scale spherical 
harmonic models of the magnetic field observed at and near Earth's 
surface. They include the 120 Gauss coefficients through degree and 
order 10 at 5-year intervals from 1945 to 1980. Linear interpolation 
gives the recommended broad-scale geomagnetic field model between DGRF 
epochs (PEDDIE, 1982); however, the work of ALLDREDGE (1987) and of 
LANGEL et al. (1988) indicates that such interpolation need not, and 
indeed should not, be taken too seriously. 
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The DGRF model for epoch 1980 consists of the Gauss coefficients of 
the model GSFC (12/83) derived by LANGEL & ESTES (1985) rounded to the 
nearest nT and truncated to degree and order 10. 

PEDDIE (1982) describes the derivation of the DGRF models for 
epochs 1965, 1970, and 1975. Each of these models is a weighted average 
of three candidate models derived by different groups of geomagnet ists : 
the GSFC (9/80) model derived by LANGEL et al . (1982); the IGS models 
derived by BARRACLOUGH et al. (1982); and the USGS models derived by 
PEDDIE & FABIANO ( 1982) . The weighted averages are rounded to the 
nearest nT and truncated to degree and order 10. The t ime -dependent 
weights assigned to the GSFC, IGS, and USGS models are, respectively, 
(0.2, 0.6, 0.2) for 1965; (0.25, 0.50, 0.25) for 1970; and (0.3, 0.4, 
0.3) for 1975 (PEDDIE, 1982). 

The DGRF models for epochs 1945, 1950, 1955 and 1960 were derived 
by LANGEL et al . (1988). The Gauss coefficients published in that paper 
are rounded to the nearest 0.1 nT, are accompanied by simultaneously 
derived SV coefficients, and include uncertainty estimates. The latter 
estimates depend upon prior assumptions including: (1) the spatial power 
spectra of the core and crustal fields used to form correlated data 
weight matrices and (2) the requirement that in regions of poor data 
coverage the models do not stray too far from a MAGSAT model projected 
back in time (LANGEL et al . , 1988). 

Spherical harmonic models of the geomagnetic field have been 
derived for about 150 years; they represent a vast number of 
observations quite well. The 8 DGRF models can be (and with general 
weights are) viewed as 8 points in a 12 1 -dimens ional sub-space: 120 
dimensions for the coefficients and one for time. The covariance 
matrices derived for the DGRF models (see section 3,1) describe 
estimated error ellipsoids around each such point (excluding 1980). My 
strategy for estimating a single steady flow is to take one of these 
points (e.g., 1980) as an initial condition and seek a trajectory whose 
projection onto the 121-dimensional sub-space passes close to the other 
7 points. The trajectory is the predicted magnetic field determined by 
upwardly continuing the results of steady, frozen-flux motional 
induction at the CMB through a SFM. If a trajectory is found whose 
projection passes close enough to (e.g., intercepts) all 7 error 
ellipsoids, then the steady flow hypothesis and the SFM/FFC earth model 
are compatible with the DGRF models and any geomagnetic data they 
adequately represent. Indeed, the steady flow hypothesis would then 
seem to be in the same class as the candidate DGRF models. If no such 
trajectory is found, then the steady flow hypothesis or the earth model 
should be suspected of significant error. If no such trajectory exists, 
then one or more of these suppositions is too poor to provide a 
quantitatively adequate explication of the indicated secular change. 

Application of the method derived in IB (or similar methods) to a 
particular set of Gauss coefficients (or original data) presents 
particular technical problems which must be solved either before or 
during the calculations. The rest of this section describes how these 
problems have been solved for the case of the DGRFs . 

2.1 Initial Conditions 

Solutions of the forward steady motional induction problem at the 
top of a FFC can exhibit an extreme sensitivity to initial conditions 
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common to many problems in deterministic chaos (IA, section 3). Initial 
conditions should therefore be selected with care. The 1980 DGRF model 
represents global geomagnetic data gathered by NASA ' s Magsat satellite 
as well as other types of data, so it was chosen to constrain the 
initial conditions by selecting t Q = 1980 when fitting either the entire 
35-year sequence of DGRF models <tf = 1945) or any sub-interval 
including epoch 1980. This is a relatively ambitious initial condition 
because 1980 is an endpoint rather than an interior point of the 35-year 
time interval, and because errors in the initial condition are expected 
to grow exponentially with elapsed time It - t Q | (see equation 
(IA.14b)). This choice will, however, enjoy an advantage when extending 
the steady motions hypothesis to make and test predictions for 
subsequent epochs . 

A complete initial condition is required. Though complete in the 
spatial domain, the DGRFs are truncated spherical harmonic models and 
are thus incomplete in the Four ier-Legendre spectral domain. The 1980 
DGRF model contains no information about Gauss coefficients of degrees 
greater than 10, so some supposition is needed to complete this model. 
I have supposed that all Gauss coefficients of degrees greater than 10 
are zero at 1980-the expectation value indicated by any probability 
distribution for these Gauss coefficients which enjoys zero mean. This 
supposition, along with the SFM/FFC model, limits my ability to falsify 
the steady flow hypothesis: if this initial condition does not allow 
derivation of a satisfactory fit to the other 7 DGRF models, there may 
be a different initial condition with a non-trivial narrow-scale core 
field which does. Because the 1980 DGRF model is a truncated version of 
model GSFC (12/83), which includes main field coefficients through 
degree and order Np = 13, coefficients of degrees 11, 12, and perhaps 13 
can be used to test the sensitivity of misfit, estimated velocity 
fields, and geomagnetic field forecasts to changes in the narrow-scale 
part of the initial magnetic conditions. Indeed, a test case showed 
that including non-trivial initial coefficients through degree 12 gave a 
very slight improvement in fit to the other 7 DGRFs. Effects of still 
higher degree initial coefficients could be investigated via Monte Carlo 
inversions of deterministic chaos-a seemingly unenlightening pursuit. 
Reliance upon non-trivial, narrow-scale initial geomagnetic conditions 
to provide a satisfactory fit may be possible, but it seems inelegant; 
indeed, contrivance of such conditions might well be dismissed as a 
fluke. Of course, with a much more realistic earth model, such design 
could reveal fine core field structure which is otherwise masked by the 
crustal field. 

2.2 Magnetic Truncation 

The spherical harmonic expansions for the input radial field (la) 
and the predicted radial field (lb) involve summations over repeated 
subscript i which ought to run from one to infinity (or to the highest 
wavenumber describing the macroscopic geomagnetic field). The same is 
true of the summations required to evaluate the square-weighted residual 
(5b) and thus to estimate the (similarly infinite dimensional) vector of 
flow parameters These sums must be truncated in practice. 

The DGRF models do not specify the Gauss coefficients of degrees 
greater than 10 (n > Ng = 10) and thus do not specify the radial field 
coefficients g^ with i > Ng(Ng+2) = 120. This does not imply, nor is it 
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tantamount to assuming, that such high-degree coefficients are zero! It 
means only that such high-degree coefficients are not known well enough, 
are hopefully not important enough, and are perhaps too influenced by 
the crustal field, to include in the DGRF models. Therefore, if the 
spherical harmonic representation of a predicted magnetic field model 
contains non-zero radial field coefficients y ^ with i > 120, such 
coefficients must not be directly included in the squa r e - we i gh t ed 
residual. This requirement can be met by truncating both g^ and y i in 
(5b) to length i rnax = 120 when evaluating and A r ^ . Such truncation is 
consistent with the assumption that high-degree DGRF cofficients are 
unknown. The possibility of using prior beliefs or biases about the 
evolution of the narrow-scale core field to constrain predicted high- 
degree coefficients is not pursued here. Such coefficients can be 
derived by extending spherical harmonic analysis of the predicted field 
past the 10 th degree and they may then be examined in detail. This was 
not done routinely; however, both the rms predicted radial field 
component, < [ B r p ( b , t ) ] ^ ^ ^ ^ , and the rms time rate of change of the 
predicted radial field, < [3 t B r p (b, t ) ] 2 > 1 / 2 , are routinely calculated to 
verify that implausible contributions from narrow-scale radial fields at 
the CMB are in fact not predicted. The former does not exceed a few 
gauss and the latter typically does not exceed several thousand nT/year. 

An alternate approach to the magnetic truncation problem is through 
the radial field weight function W(a,t) and weight matrices W(a,t) (or 
the general weight matrices Q(a,t)). The objective is to match the 
evolution of the broad-scale part of the radial field alone, so the 
weight function must reflect uncertainties in the broad-scale radial 
field represented by low-degree coefficients alone. Then (possibly 
unknown) radial field coefficients of degrees greater than some maximum 
degree (here 10) must be assigned zero weight. To do so one may set the 
W — with i or j (or both) greater than 120 equal to zero. Similarly, if 
the objective is to match the evolution of the broad-scale scalar 
potential, then direct contributions to the square-weighted residual 
from high-degree coefficients must be assigned zero weight. This can be 
accomplished by truncating each Q(t) to a 120x120 matrix. Truncation 
of W (or fl) has the same effect” on the square-weighted residual as 
truncation of both g and 7 . Both approaches avoid including the direct 
effect of non-trivial predicted high-degree coefficients in the square- 
weighted residual. 

Truncation of weight matrices does not necessarily justify 
truncation of covariance matrices used to derive them. Ideally, 
complete covariance matrices would be inverted for complete weight 
matrices. Truncation of the complete weight matrices would then target 
the broad-scale objective. 

In practice, all g^ above degree 10 are taken to be zero only to 
start the iteration scheme (j = 0 in equation (7)) . They are also taken 
to be zero at time t Q when an initial condition is needed to prime the 
forward motional induction scheme. The latter is a finite difference 
representation of equation (0) which induces SV by steady numerical 
advection of the radial field on a 2 °x 2 ° regular mesh using 10 ^ year 
time steps. (The adequacy of a 2 mesh has-been confirmed by comparing 
test results with those obtained using a 1 mesh; the accuracy of the 
first-order time-stepping has been confirmed using a second-order 
scheme) . The resulting numerical solutions of the forward steady 
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motional induction problem implicitly represent the predicted magnetic 
field through very high degree: the azimuthal Nyquist wavenumber is 90 
on the 2°x2° mesh. 

This fine resolution allows key effects of the predicted narrow- 
scale core field on both the solution of the forward problem and the 
iterative solution of the inverse problem to be retained without having 
to calculate predicted coefficients above degree 10. Such high-degree 
coefficients would represent narrow-scale core field structure induced 
by the interaction of the steady fluid motion with the radial magnetic 
flux density (VOORHIES, 1984; 1986a; 1986b) . Because a steady fluid 
flow also interacts with the narrow-scale field to induce secular change 
on all spatial scales, the predicted narrow-scale fields eventually 
influence the predicted broad-scale field. Specifically, predicted 
high-degree coefficients eventually influence predicted low-degree 
coefficients (see untruncated equations IB.15-16b). Because predicted 
low-degree coefficients determine the square-weighted residual, this 
indirect effect of predicted high-degree coefficients must be included 
to accurately measure misfit. 

Only the y ^ through degree 10 contribute directly to the square- 
weighted residual, so only the through degree 10 contribute 
directly to (5b) (see (IB. 9c , 11c and 19)). It follows that linearized 
attempts to minimize the square-weighted residual need only account for 
the first 120 rows of A (VOORHIES, 1986b) . This is also true for 
attempts to solve the non-linear inverse problem by iterative solution 
of the linear problem. However, as noted above, the indirect effects of 
predicted high-degree coefficients on the predicted low-degree 
coefficients must be retained when solving the non-linear problem. This 
requirement can be fairly well, if not strictly, met by several 
techniques-each of which demands very accurate calculation of the total 
predicted radial field component needed in (4a-c). 

The physical technique adopted tracks B (b,t) in the physical 
domain via finite-difference solution of equation (0) on a fine mesh. 
Values of B r p(b,t) at the mesh points are periodically stored for use 
with (4b). Only the low-degree coefficients needed to evaluate the 
residuals and 53^g(j) are extracted from B r p* The surface integrals 
giving P and Q (IB. 17) are evaluated numerically, so both broad and 
narrow scales in the predicted field contribute to the new A(j) needed 
for the next deep iteration (7). One pseudo-spectral technique would 
solve (0) in the physical domain, but periodic spheiical harmonic 
analysis of B r p(b,t) for the r^(t) would extend through very high degree 
and order (say 90). These would then be used in (4a) with X — ^ and 
Y i_jk evaluated through extremely high wavenumber k (say k < 8,280), but 
only i < 120. A purely spectral technique based on (4a) would track 
evolution of the r^(t) through very high degree using (3), the X — ^ and 

from (IB.16a-b) evaluated to extremely high wavenumbers i and k, 
ana, say, T^t+St) * T^(t) + (t ) ]6t with very small 5t. Efficient 
evaluation of the X — ^ and Y— ^ has been achieved by K. A. Whaler and D. 
Winch (1987, personal commun icat ion ) . Spectral methods may be more 
efficient computationally, but I prefer the simple, direct physical 
method (which might be readily adapted to an aspherical CMB) . 
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2.3 Velocity Truncation and Damping 

The original method of VOORHIES (1984; 1986a) yields steady flows 
which fit SV models at the CMB under the assumption that SV coefficients 
above degree are zero. It was found that the st r eamf unc t ion should 
include coefficients through degree 12 or more in the case of purely 
toroidal flow. For combined flows, the preferred solution (C 4) included 
both streamf unct ion and velocity potential coefficients through N T = 

= N v = 12. The second method (VOORHIES, 1986b) yields steady flows 
which fit SV at Earth's surface without assuming zero high-degree SV; 
however, omission of non-linear aspects of the inverse problem and the 
avoidance of damping prohibited meaningful estimation of flow 
coefficients above degree 5. Estimates with N v > 5 gave smaller linear 
residuals, but larger genuine, non-linear, residuals obtained via 
forward solution of the steady motional induction problem, spherical 
harmonic analysis of B r p(b,t), upward cont inutat ion , and comparison with 
the DGRFs . 

After some experimentation with new algorithms, I chose to truncate 
the velocity field parameterization at degree and order N v = 16 and to 
use non-zero damping A^. Then 288 and 288 (3^, or 576 velocity field 
parameters are determined; these parameters are not entirely free 
when either A^ or Ag (or both) are non-trivial. Such truncation is 
satisfactory when equation (7) is solved with large to moderate A^ and 
J^ 0 = 0. The latter selections also yield solutions i; ( j ) which converge 
in a few iterations. When A^ is too small, shallow iteration using the 
initial A(j = 0) matrices did not appear to yield convergent solutions, 
indicating a failure of the linearization in the case of vigorous, 
complicated flow. Moreover, deep iteration with small A^ sometimes gave 
£{j) converging towards local minima characterized by greater square- 
weighted residuals than solutions of lesser roughness. This is to be 
expected when non-linear aspects of the inverse problem are not 
overwhelmed by a fixed prior bias ^ = 0 (towards zero solution norm) . 
Physically, the spectral mapping of a fixed prior bias towards zero norm 
can inhibit discovery of solutions with (a priori, improbably) tall 
peaks in the surface kinetic energy density (SKED) spectrum (VOORHIES, 
1984, 1986a) that offer a smaller misfit at equal or lesser norm. 
Linear stability calculations (CHANDRASEKHAR, 1961) lead me to suggest 
that tall peaks in the SKED spectrum are in fact quite plausible. 

The learning algorithm does not suffer so much from this drawback 
because it allows iterative growth of tall spectral peaks which may be 
needed to reduce misfit. It was used to derive rougher flows and 
tighten the upper bounds on the apparent error in the earth model. 
Examination of SKED spectra indicates N v = 16 is acceptable for flows of 
moderate roughness derived using the learning algorithm. The learning 
algorithm need not produce convergent £ ( j ) (unless A^(j+1) is set to 
infinity), so increasing N v could reduce the misfit for the roughest 
flows explored. For very large N v , the learning algorithm yields flows 
with SKED spectra which must in principle still fall off faster than n”^ 
at high n; however, the amplitude of the spectra will depend on both the 
convergence factors (A^(j)) used in, and the solutions obtained by, all 
previous iterations. Such hysteresis effects may occur in other 
applications of the learning algorithm. 
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2.4 Integrations 

The double time integrals in equation (7) can be approximated 
numerically in many ways. I chose a simple method based on accepting 
the suggested linear interpolation between DGRF models at face value. 

The interval t Q to t f (e.g., 1945 to 1980) is subdivided into L = 7 
sub- intervals . For each sub-interval 1 , let t^ denote the mid-point, 
t^_ the beginning, t^ + the end, and let At^ = t 1+ - t^_ denote the 
duration (5 years). The time rate of change of the radial field is 
approximated by a temporally constant, but spatially variable, SV in 
each such sub-interval in accord with the suggested linear 
interpolation. The t ime -dependent A matrix is the upward continuation 
of ‘the concatenated P and Q matr ices - ( 4c ) . Elements of P and Q (IB. 17a- 
b) are linear in the~radial field at the CMB, so the integral of any 
such element over any sub-interval At^ is approximated by the value of 
that matrix element at the sub-interval midpoint t^ multiplied by At ^ . 
These values were estimated by substituting the arithmetic mean of 
B r p(b,t) at times fc i+ and t^_ into equations (IB.17a-b) and evaluating 
the integrals over the CMB numerically using a simple double Riemann sum 
approximation on a 2°x2° mesh. 

The dummy time integral of the A matrix is thus approximated by its 
Riemann sum 


1 + 


A (X ) ( t ) = j A { X ) dx « A ( 1 1+ ) = J A ( X ) dx - X A(t k )At 


k=l - 


:8a) 


and the elements of 83 x g at time t are approximated by 

t t 1+ 

8d T g(T) (t) a I 6a T g(x)dx « 83 T g(X) (t 1+ ) 83 x g(x)dx 


« X 53 x g(t k )At k 
k=l 

where 


(8b) 


^t^i^k^ = tl9i(tk + ) ST i ( t k— ) 1 1 H At k ) 


-1 


(8c) 


If the initial condition were at 1945, then k would be incremented from 
1 (sub-interval 1945 to 1950) to 1 (sub-interval t^_ to t^ + ). With the 
initial condition at 1980 (the end of the interval), k is actually 
decremented from 7 ( sub- interval 1975 to 1980) to 1. However, it is 
convenient (and thermodynamically permissible) to imagine time running 
backwards in this dissipationless (and thus isentropic) problem. Then 
relabeling the sub-intervals still allows (8) to be used. 

Because A and 83^g are treated as steady during each sub- interval 
in the dummy time integrals, it is convenient to further approximate the 
weight matrix elements W • ■ as constants during each sub-interval. Then 

J- J m 

the double time integral of A WA in equations (7) is approximated by 
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(9a) 


L 1 1 

A ( X ) T W ( t ) A (T) - I {[ I A(t k )At k ] W(t x ) [ Z A(t k )At k ] }At : 
1=1 k=l = = k=l = 


and the double time integral of A^Wfx^g is approximated by 


L 1 1 

A(x) T W(t)6a x g(x) - Z {[ Z A(t k )At k ] w(t x ) [ Z 5a x g(t k )At k ] JAt,. (9b) 

1=1 k=l = = k=l 

The finite time element approximations (9a-b) seem particularly 
appropriate for initial iterations (j =0) using the DGRFs-provided the 
weight matrices are roughly constant. If they are not, it would be 
better to use a more precise technique (e.g., select L > 7). This has 
not been attempted because the best method for interpolating the DGRF 
models and their weight matrices is not known. (A more accurate method 
for evaluating the dummy time integral for j > 0 which exploits the 
results of the forward solution has recently been adopted) . 

No allowance for the polar caps was made in the spatial numerical 
integrations used to calculate elements of P and Q, hence A. These 
matrix elements are thus accurate to only about two parts in 10 4 . 
Indeed, downward continuation of the 1980 DGRF to the CMB, evaluation of 
B r (b,1980) on the mesh, spherical harmonic analysis through degree 10 
via numerical integration, and subsequent upward continuation yielded 
rms errors in the radial field of 6.77 nT (about 2 parts in 10 4 ). This 
error is due primarily to inexact reproduction of zonal (i = n 2 ) 
coefficients and to the large axial dipole in particular. Although this 
error is smaller than the minimum 25 nT rms uncertainty in the pre-1980 
DGRF radial field models (see section 3.2) , a polar correction was 
included in the spherical harmonic analysis of the predicted field for 
the Y i at all DGRF epochs so as to calculate the misfit precisely. This 
reduced the error from 6.77 nT to 0.11 nT in the test case. 

Small errors resulting from imperfectly exact integrations are 
unimportant compared with expected errors in the underlying physical 
assumptions. For example, it is doubtful that the radius of the core is 
known to better than 1 part in 10^: one survey of seismic estimates 
indicated the core radius to be 3485 ± 5 km (VOORHIES, 1984). Because 
errors in the core radius 8b (and larger errors associated with the 
omission of core ellipticity) amplify as [ a / ( b+8 b ) ] n+ 1 when B r is 
extrapolated from Earth's surface to the CMB, they should exceed the 
spatial integration errors. Fine tuning of the spatial integration is 
not too difficult, and may be needed when core asphericity is included. 
(Indeed, extremely accurate spatial integrations based on a quartrature 
rule have recently been developed and implemented) . 

Tests conducted with synthetic field models show that the numerical 
methods do not contribute significantly to the misfit, A r 2 . These tests 
also reveal that recovery of steady flows used to synthesize test case 
secular change models, and thus the effective elimination of misfit, 
requires extremely small damping X Nevertheless, I can but establish 
upper bounds on the error in the earth model because of the uncertainty 
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and 


in narrow-scale initial magnetic conditions, the finitude of N v , 
limitations imposed by either fixed bias or learning algorithms, 

3. DERIVATION AND IMPLICATIONS OF WEIGHT MATRICES FOR THE DGRF MODELS 

3.1 Derivation of the Weight Matrices 

Derivation of weight matrices for the 1945, 1950, 1955 and 1960 

DGRF models proved straightforward. The total information matrices for 
these models were generously provided by R. Langel and T. Sabaka (LANGEL 
et al., 1988). These 200x200 matrices were inverted to obtain the full 
covariance matrices for these DGRF models. The latter were stripped of 
main field-secular variation (MF-SV) and SV-SV coefficient covariances 
because the main field weighting scheme was selected. The square roots 
of the diagonal elements of the resulting 120x120 main field covariance 
matrices, V{t), are the uncertainty estimates for the main field 
coefficient's published by LANGEL et al . (1988). The elements V^tt) 

were multiplied by [ n ( i ) + 1 ] [ n ( j ) + 1 ] to obtain the epoch-dependent 
elements of the covariance matrix for the input radial field 
coefficients, E i j(t n ). These matrices were used to calculate the 
squared uncertainty estimates for the DGRF models at the mesh points of 
the 2°x2° regular grid using the inverse of equation (IB. 10). The 
resulting OB (a,t n ) f the square root of the inverse weight functions 
W(a,t)~ 1//2 , are the same as the uncertainty estimates contoured by 
LANGEL et al. (1988, figures 2c, 3c, 4c, and 5c). As expected, the 

weight function is large (heavy) in areas of dense data coverage (e.g., 
Europe) and small (light) where data are sparse (e.g., the SE Pacific 
Ocean and portions of the Indian Ocean) . The weight matrix elements W i • 
at DGRF epochs 1945, 1950, 1955, and 1960 were evaluated by numerical 

integration of (IB. 12a). 

A very different method was used to derive weight matrices for the 
DGRF models at epochs 1965, 1970, and 1975. Recall that a weighted 

average of the three candidate models defines the DGRF models at these 
epochs. The weight factors are supposed to reflect the apparent 
accuracy of the candidate models (PEDDIE, 1982) . This weighted average 
was taken as defining the expected value of the Gauss coefficients and 
thus the expected value of functions thereof. 

Let the vector of radial field coefficients for candidate model k 
at epoch t n be denoted by its elements g.j^(t n ) . Then the radial field 
coefficients for the DGRF model at epoch t n are 

g i ( t n ) - 9i (t n> s w ij k(t n ) 9j k(t n ) (10) 

where the approximation represents roundoff error. The weight tensor 
v .^(t n ) depends on epoch t n , candidate model k and the coefficient of 
interest i; its projection at a given t n and fixed k is a matrix which 
is diagonal in the coefficient subscripts i and j. At epochs 1970 and 
1975 all candidate models are complete through at least degree 10; w— 
is then the non-zero constant w^ for i = j < 120 and zero for all other 
( i t j ) . Higher degree coefficients, including the GSFC (9/80) 
coefficients of degrees 11 through 13, are assigned zero weight; this 
does not imply that they are zero. 

The 1965 IGS candidate model includes coefficients only through 
degree 8 ( i < 80) . Inspection reveals the 1965 DGRF coefficients of 
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degree 9 and 10 (81 < i < 120) to be the rounded, equally weighted 

averages of the GSFC and USGS candidate models. This is in accord with 
the assignment of zero weight to the undetermined degree 9 and 10 
coefficients of the 1965 IGS model. These coefficients are not assumed 
to be zero here, so the 1965 weight tensor elements depend upon 
coefficient degree. Note that the assignment of zero weight to certain 
geomagnetic coefficients is a geophysical judgment; it is not a 
mathematical result obtained by inverting infinite dimensional 
information matrices with elements either determined or bounded by 
either data or prior information. 

Equation (10) is here taken as defining the expected value of the 
DGRF model coefficients. This expected value is a linear combination of 
the candidate model coefficients: the weighted average with weights 

inversely related to the uncertainty ascribed to the candidate models. 
Linear combinations of DGRF coefficients are thus linear combinations of 
the candidate models, so the expected value of the error of any linear 
function of DGRF coefficients is zero (apart from roundoff error) . For 
example, the expected broad-scale radial field component 

E{B r (a,t n )} = E{g i k (a,t n )S i (0,()»} = [w i:j k (t n ) g j k (t n )]S i 

- g i (t n )s 1 - g i (t n )s i (lla) 

is that given by the DGRF coefficients (to within roundoff error) . The 
individual coefficient deviations 


d i k (t n ) * 9i k ( t n ) - 9i ( V (llt>) 

have been calculated using the published candidate and DGRF model 
coefficients. The expected error in the broad-scale radial field 

E{5B r (a,t n )} = E{d i k (a,t n )S i (0,<t)) } = [w i:j k (t n ) d j k (t n )]S i 

* [«i<t n > - g i (t n )is i - o uio 

is indeed zero to within the expected cumulative (rss) roundoff error of 
0.7 nT per Gauss coefficient. Note that non-linear functions of DGRF 
coefficients like the scalar geomagnetic intensity or the absolute flux 
linking the core (BENTON Sc VOORHIES, 1987) are non-linear combinations 
of the candidate models and need not have zero expected error. 

The expected value of the squared error in the broad-scale radial 
field is needed to calculate the weight function (IB. 10). This is 

E{ [8B r (a,t n ) ] 2 } = E{ [S i d i k ] [dj k Sj ] } = S i E{[d i k ] [dj k ] } Sj 

= Si E i:j (t n ) Sj = (OB r (a,t n )] 2 = W(a,t n ) _1 . (12a) 

Cross correlations between candidate models (different superscripts) and 
different epochs (t n ) are omitted in accord with (10). ^For each t n and 
each k, introduce the diagonal matrix square root of Wij K (t n ), CO: 


®im (k '’ t n ) “mj (k;t n ) 


W 


11 


•(t n > 


(12b) 
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with no sum over k or n . The expected value of the squared uncertainty 
in B r , [OB r (a,t) ] 2 , is 

S i E ij s j = s i ( [Wii k d : k ] [co jm k d^]} Sj = Wfa.t)" 1 (12c) 

and the covariance matrix for the radial field coefficients is 


Eij = [L ± CO,,* df] [1^ to jm k dj*] = I k Eij 1 


(12d) 


where the sum on k is performed after the relevant i by j matrices have 
been formed. The covariance matrix is thus a weighted sum of dyadics . 


to be weighted in 
the elements of the 
squared uncertainty 
elements (12d) have 


The © matrix allows the coefficient deviations 
accord with the DGRF averaging procedure. Because 
weight tensor w are normalized to unity, both the 
in the radial ?ield (12c) and the covariance matrix 
the correct units of nT 2 . 

The radial field covariance matrix elements 
using the published Gauss coefficients and 
uncertainty (conditional standard deviation) of 
coefficients is the ratio of the square root of the diagonal elements of 
E— to n(i)+l. At epochs 1970 or 1975 the elements of symmetric E are 

This is not true at epoch 1965, when the 


E i j 
( 1 2d ) . 


were calculated 
The estimated 
the DGRF Gauss 


"ID 

of the form 


Eij = i k di^Ai, 

Efj with both i > 90 and j < 80 reflect both the non-zero probability 
(1 - (0.5x0. 2) = 68.4%) that they are unknown and the naive 
expectation that unknown elements of E— are zero. Other expectations 
could well overestimate the uncertainty in the DGRF models. Such poorly 
determined cross -correlat ions between different coefficient deviations 
are assigned reduced weight in my attempts to weight these evolving 
broad-scale geomagnetic field models. 


The squared uncertainty estimates W(a,t n ) ^ were evaluated at the 
mesh points of the regular 2°x2° grid using (12c) (the inverse of 
(IB. 10)) . This quantity is positive at all such points; however, at a 
few (4 or 5) points it is less than 1 nT 2 , indicating possibly 
fortuitous outstanding agreement between candidate models. This can be 
understood in terms of contours of perfect agreement between candidates: 
contours where the first and second models agree can intersect contours 
where the second and third models agree. The points of intersection 
enjoy singular weight. More typically, the weight function again proved 
to be heavy in regions of dense data coverage (e.g., Europe) where the 
candidate models should agree, and light in regions of sparse data 
coverage (e.g., the SE Pacific Ocean and portions of the Indian Ocean) 
where the candidate models can differ more easily. The weight matrix 
elements W— at DGRF epochs 1965, 1970, and 1975 were then evaluated by 
numerical integration of (IB. 12a). 

The weight matrices used in equations (9) and thus the solutions of 
(7) must reflect the typical uncertainty in the broad-scale radial field 
during the sub- intervals between DGRF epochs. These were estimated by 
numerical integration of equation (IB. 12a) with the weight function 
taken to be the inverse of the arithmetic mean of the DGRF inverse 
weight functions at the sub-interval endpoints. Between DGRF epochs the 
squared uncertainty in the radial field [<JB r (a,t)] 2 is thus treated as 
if it were the (spatially variable) constant equal to the average of its 
values at the sub-interval endpoints. For the interval 1980-1975, this 
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procedure was used with the additional assumption that the uncertainty 
in the initial radial field at epoch 1980 is negligible. The estimated 
rms uncertainty in B r (a,t) at 1980 for the GSFC 12/83 through degree and 
order 10 is much smaller than that for the 1975 DGRF (16 nT vs. 58.30 
nT) , so this is a fair approximation. 

For general weights, the Cl ij(t n ) are use d in place of the W— (t n ) 
described above. The derivation of Cl matrices and the effect of using 
them are described in the Appendix (also see the Appendix of Part IB) . 

3.2 Weighted Variance in the 1980-1945 DGRF Models 

There are many quantities of geophysical interest which can be 
determined from the DGRF models and their covariance matrices. Of 
immediate interest are the more easily interpreted scalar quantities 
which describe the signal and the noise in these models as viewed from 
the standpoint of one attempting to model geomagnetic secular change 
relative to the 1980 MAGSAT epoch. 

One such quantity is the root -mean - square change in the radial 
magnetic flux density relative to 1980 averaged over Earth's surface 

AB r (a,t) = <[B r (a,t> - B r (a, 1980 ) ] 2 > 1/2 (13) 

(VOORHIES, 1986b). Another is the rms uncertainty in the radial field 
obtained from the DGRF covariance matrices described in section 3.1 


OB r (a,t) = <[aB r (a,t) ] 2 > 1/2 . 

(14a) 

Such integral measures have many representations; 
square of (14a) is 

for example, the 

[OB r ] 2 = <W(a, t) _1 > 

(14b) 

= <s k E kl s l > 

(14c) 

120 


= I {E ii /[2n(i)+l] ) 

( 14d) 


i = l 
120 

= I { ( og ± ) 2 / [2n(i)+l] ) ( 14e) 

i = l 

10 (n+1 ) 2 n 

= I I { [Og n m ] 2 + [Oh n m ] 2 } (14f) 

n=l 2n+ 1 m-0 

where (14d) follows from (14c) by orthogonality of the spherical 
harmonics, the og^ in (14e) are the estimated uncertainties in the 
radial field coefficients, and (Cg n m , O n m ) in ( 1 4 f ) are the estimated 
uncertainties in the Gauss coefficients of standard notation. The ratio 
AB r (a, t ) /OB r (a, t ) provides a crude dimensionless measure of how 
significant a change in the radial field is called for by the DGRF model 
at time t relative to the 1980 DGRF model. 
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A better measure of the secular change signal in a DGRF radial 
field model relative to the 1980 DGRF is the weighted version of (13) 

o 7 1/2 

5 0 (a,t;t o ) = <[B r (a,t) - B r (a, 1980) ] 2 W(a,t)> (15) 

1/2 

= ((Cg i (t k ) - gi ( 1980 ) ] W i j(t k ) [g-j(t k ) - gj (1980) ] }/47t) 

with sums over i and j running from 1 to 120. This quantity measures 
the secular change signal at time t (relative to 1980) in units of 
estimated uncertainty (noise). It is the instantaneous rms signal-to- 
noise ratio. The 4 it semi-normalization allows direct comparison of 8 q 
to AB r /OB r . The quantity 4 tc8q^ is the instantaneous weighted variance 
in the secular change relative to 1980. 

Time integration of the instantaneous weighted variance yields the 
cumulative weighted variance 

2 t 2 

41C[A 0 (a, t;t Q ) ] = An j [8Q(a f t;t 0 )] dt 

fc o 

t 2 

= 47tJ <[B r (a,t) - B r (a, 1980) ] W(a,t)> dt . (16) 

fc o 

At t = t £ this is the total weighted variance in the secular change of 
the radial field targeted for reduction. This is the square-weighted 
residual resulting from the absurd hypothesis of no secular change: a 
constant broad-scale radial geomagnetic component equal to that at 1980. 

The method of interpolating both DGRF models and the weight 
function should be specified to evaluate the time integral in (16). The 
recommended interpolation is linear, but the working hypothesis 
indicates interpolation via steady motional induction at the top of a 
FFC bounded by a SFM. Because linear interpolation is not taken too 
seriously, and because DGRF models are provided only at discrete epochs, 
the cumulative weighted variance in the secular change of the radial 
field called for by the DGRF models is taken to be 

n 

4n[A 0 ' (a, t n ;t Q ) ] 2 = An I [8 0 ' (a, t k ; t D ) ] 2 (17) 

k=l 

n 

= I [g i (t k ) - (1980 ) ] W ij (t k ) [g-j(t k ) - g j ( 1980 ) ] 

k=l 

where t n and t k are DGRF epochs. The cumulative weighted variance (17) 
is conveniently independent of the interpolation method and is 
dimensionless, as is the total weighted variance An [ Aq ' ( a, t n =tf ; t Q ) ] 2 . 

When the general weights based on Q — (t R ) instead of w ij(t n ) are 
used, the cumulative weighted variance is 

n 

/[A q * ( a, t n ;t Q ) ] 2 = / L [8 0 * (a, t k ;t D ) ] 2 (18) 

k=l 
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n 


= I [g ± < t k ) - g ± (1980) ] fl*ij(t k ) [g-j (t k ) - g j ( 1980 ) ] 


where the semi -normalization factor / is taken to be 120 (see Appendix) . 
This choice allows direct comparison between 5 q and 5 q' and between Aq * 
and Aq ' . 

Table 1 lists, as a function of t n , the rms uncertainty in B r , O B r 
(14); the rms change in B r , AB r (13); the ratio AB r /OB r ; 5 q' from (15); 
and the semi -normal ized cumulative weighted variance, A n ' ^ ( 17 ) . Table 
2 lists a B r ; AB r ; the ratio AB r /o B r ; 8 q ; and Aq z appropriate to 
general weights. Note (a* B r ) 2 = E ii / [ 2n ( i ) +1 ] = (CJB r ) 2 + (11.24 nT) 2 
for epochs 1965, 1970, and 1975 (see Appendix) . 

Table 1 shows that a B r is about 51 nT for the 3 averaged DGRF 
models. Though little or no satellite data were available at epoch 
1975, the decline of AB r going back from 1975 to 1965 likely reflects 
the heavier weight ascribed to the IGS candidate model at the earlier 
epochs. Note oB r is but 26 nT at 19 6 0-poss ibly because of improved 
methods used to derive the 1960 DGRF and its covariance. The rms change 
in the expected broad-scale radial field AB r between 1945 and 1980 is 
over 1900 nT . The typical rate of change of B r was about 55 nT/yr, but 
was more rapid later in the in terval-about 64 nT/yr between 1975 and 
1980-in part because of the accelerating decline of the axial dipole. 
The ratio AB r /oB r is always less than Sq', so the weights increase the 
signal magnitude (by an rms factor of 1.9). It follows that regions of 
heavier weight (higher accuracy) correlate with regions of faster 
secular change (larger AB r ). This might be explained by the relatively 
low areal data density over the Pacific Ocean where secular change is 
thought to be slow. 

The changes in the rms signal magnitude 8 q' as one looks back in 
time reflect competition between increased change in the radial field 
due to secular change (increased unweighted signal AB r ) and decreased 
weight (increased noise level GB r ) prior to 1960. The latter is likely 
due to a sparser data distribution before the International Geophysical 
Year and the advent of satellite geomagnetomet ry . Caution in accepting 
this explanation is advised because the prior information used in 
deriving the 1945-1960 DGRF might reduce the apparent unweighted signal 
magnitude in regions of poor data coverage. The semi -normalized total 
weighted variance is 12,217.8; the square root of this number, 110.5, is 
the total rms signal-to-noise ratio. 

Table 2 shows that the general weights yield much larger values of 
the rms signal-to-noise ratio and the semi -norma 1 i zed and cumulative 
weighted variance than do the radial field weights. This is not too 
surprising because the different weights correspond to two different 
objective functions (see Appendix of Part IB). Indeed, radial field 
weights should indirectly filter out the effects of D, X, Y, and H 
geomagnetic data which help determine the scalar geomagnetic potential. 
The semi -norma 1 i zed total general weighted variance is 57,548.9, 
indicating a total general rms signal-to-noise ratio of 239.89. 

The semi -normalizations allow comparison of weights; they are not 
intended to be misleading. With t = 1980 and tf = 1945, the non- 
normalized total weighted variances are 4 tc [ Aq ' (a, t f ; t Q ) ] 2 = 153,533 and 
120 [Aq (a,tf;t 0 )]^ = 6,905,832. The ratio of these numbers is about 
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TABLE 1 


Select Properties of the DGRF Models and Uncertainty Estimates 
Based on Radial Field Weights 


fc n 

(A.D. ) 

GB r (a,t) 

(nT) 

AB r (a, t ) 

(nT) 

AB r /oB r 

V (a ' t n ;t o ) 

V< a 't n 

1980 

0.0 

0.0 

- 


- 

1975 

58.30 

318.9 

5.47 

13.64 

186.0 

1970 

47.47 

612.7 

12.91 

26.01 

862.5 

1965 

45.42 

899.2 

19.79 

51.96 

3562.8 

1960 

25.67 

1157.3 

45.08 

61.17 

7304.4 

1955 

65.87 

1427.2 

21.67 

64.35 

11445 . 1 

1950 

127.01 

1674.9 

13.19 

19.24 

11815.1 

1945 

160.89 

1929.3 

11.99 

20.07 

12217.8 


TABLE 2 

Select Properties of the DGRF Models and Uncertainty Estimates 
Based on General Weights 


*n 

(A.D. ) 

a B r (a,t) 

(nT) 

AB r (a,t) 

(nT) 

AB r /a B r 

So 

A 0 ^ a,t n' 

1980 

0.0 

0.0 

- 

- 

- 

1975 

59.37 

318.9 

5.37 

42.09 

1771.8 

1970 

48.78 

612.7 

12.56 

79.95 

8164.3 

1965 

46.79 

899.2 

19 .21 

119.17 

22366.9 

1960 

25.67 

1157.3 

45.08 

163.40 

49067 .6 

1955 

65.87 

1427.2 

21.67 

72.73 

54357 .7 

1950 

127 .01 

1674.9 

13.19 

35.69 

55631.6 

1945 

160.89 

1929.3 

11.99 

43.78 

57548.6 
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2.222%; their square roots are 391.833 and 2,627.89; and the square root 
of their ratio is about 14.91%. This latter figure suggests that the 
radial field weights filter out 85% of the geomagnetic observations. In 
this sense it takes 67 of the geomagnetic data modeled to affect the 
expected radial field in the same way as 10 of the radial field data 
modeled. This is not too surprising since Z data alone could determine 
the scalar geomagnetic potential, but are uncommon. 

3.3 The Target Residual Variance 

The definition of the cumulative weighted variance in the DGRF 
models (17) requires the cumulative square-weighted residual relative to 
the DGRF models to be redefined as 

2 n 2 
4rc[A r ' (a, t n ;t 0 ) ] = 4 tt Z [8 r ' (a , t k ; t Q ) ] 

k=l 

n 

= 2 [ffi<t k ) - Y i (t k )] T W i j (t kUg j (t k ) - Y-j <t k ) ] (19) 

k=l 

with i < 120. This measures how well the predicted radial field fits 
the DGRF models . At t n = t^ it is the total square-weighted residual . 
For t Q = 1980 and t f = 1945, if the residual at each pre-1980 epoch is 
one estimated standard deviation (10'), then the semi-normalized square- 
weighted residual 8 r ^ is (lcr')^ at each epoch and the semi -norma lized 
total square-weighted residual A r ' (a, t f ; t } ^ is 7. 

Replacing the primes with asterisks, W with Cl * , and 4 n with 120 in 
(19) yields the cumulative generalized square-weighted residual 

* 2 n 2 

120[A r *(a,t n ;t o )] = 120 I [5 r *(a, t R ; t Q ) ] 

k=l 


n 

= 2 fg i (t k ) - Y i (t k ) ] T Q i:j (t k ) [ 9j (t k ) - Yj (t k ) ] . (20) 
k — 1 

At = t f this is the total generalized square-weighted residual. For 
to = 19 *® and t f * = the residual at each pre-1980 epoch is lo*, 
then 8 r is (lo ) 2 at each epoch and A r * (a, t f ; t Q ) 2 is 7. This is the 
total square-weighted residual left by any trajectory whose projection 
onto the DGRF sub-space grazes each of the 7 error ellipsoids. Other 
trajectories can of course hit the target semi -normal i zed square- 
weighted residual of 7 provided residuals greater than lo* at some 
epochs are balanced by errors less than lo* at other epochs. 

If one accepts the uncertainty estimates derived for the DGRF 
models and the radial field weighting, then the target square-weighted 
residual is 7 out of 12,217.8. Then 99.9427% of the total weighted 
variance in the evolution of the broad-scale radial field indicated by 
the DGRF models is signal to be fit. The residual 0 . 0573 % may well be 
noise and defines the target normalized square-weighted residual. Its 
square root, 2.394%, is the typical noise-to-signal ratio m the secular 
change of the radial field described by a DGRF model relative to 1980. 
The typical signal-to-noise ratio is ( 12, 217 .8/7 ) 1/2 or 41 . 78 . Typical 
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( rms ) residuals of lo', 2 o', or 30' correspond respectively to total 
square-weighted residuals of 7, 28, or 63 out of 12,217.8. 

If one prefers the general weights, then the target total square- 
weighted residual is 7 out of 57,548.6. Then 99.9878% of the total 
generally weighted variance in the evolution of the broad-scale scalar 
geomagnetic potential indicated by the DGRF models is signal to be fit. 
The residual 0.0122% is noise and defines the target normalized square- 
weighted residual. Its square root, 1.103%, is the typical noise-to- 
signal ratio in the secular change of the scalar geomagnetic potential 
described by a DGRF model relative to 1980; the typical signal-to-noise 
ratio is 90.67. This is over twice that indicated by the radial field 
weights. Typical (rms) residuals of lo*, 2a*, or 3a* again correspond 
respectively to total generalized square-weighted residuals of 7, 28, or 
63 out of 57,548.6. 

In my opinion, residuals less than 0.674a suggest the unnecessary 
fitting of noise; residuals exceeding lo suggest a flawed hypothesis. 
The hypothesis to be tested is that the secular change required by the 
DGRF models can be explained by a source-free mantle surrounding a 
frozen-flux core in (at least temporarily) steady surficial flow. 
However, the SFM/FFC model is but an approximation, the core flow is not 
strictly steady, the initial conditions are not very well known due to 
ignorance of the narrow-scale core field, numerical approximations are 
used in applying the method derived in section 3, and my uncertainty 
estimates are but estimates. Moreover, normalized residuals of 10% 
indicate a model of secular change which is accurate to first order. 
Yet such 10% residuals indicate a typical error of 9.06a (or 4.18a r ') 
per DGRF epoch . 

In my opinion, if the hypotheses of steady flow at the top of a FFC 
surrounded by a SFM can be used to derive a fluid velocity field which, 
via steady motional induction, reproduces radial field models to within 
normalized residuals of 10%, then it has some qualitative merit. 
Residuals less than 3a' (6.52a*) indicate appreciable qualitative merit. 
Residuals less than 2a' (4.34a*) indicate some quantitative merit. 
Residuals less than la' (2.17a*) indicate appreciable quantitative 
merit. Residuals less than la* demonstrate adequate quantitative merit. 
However, any residuals greater than la* indicate a need to improve 
either the earth model (by eliminating one or more of the underlying 
assumptions), or upon the means used to test it. 

4. SUMMARY 

Application of the method developed in paper IB to the Definitive 
Geomagnetic Reference Field (DGRF) models [IAGA, 1988] has been 
described. The numerical methods used and the derivation of weight 
matrices for the DGRF models were outlined. The DGRF models are 
truncated spherical harmonic representations which contain no explicit 
information on any narrow-scale structure of the geomagnetic field. 
Attempts to fit such models should not be penalized for predicting the 
existence of such structure, provided its amplitude is not unreasonably 
large. The derived weight matrices reflect this fact, were used to 
assess the information on secular change contained in the DGRF models, 
and were further used to estimate the tolerable residuals which could be 
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produced by a quantitatively acceptable explication of definitive 
secular change. 

Preliminary results of applying the method have been presented 
elsewhere (e.g, VOORHIES, 1987b, c; 1988a, 1989b); final results are to 
be presented in part ID. 

Application of the method to the DGRF models immediately encounters 
a formidable obstacle: incomplete initial magnetic conditions in the 
Fourier -Legendre spectral domain. The Gauss coefficients above degree 
10 were therefore treated as if zero to complete the initial condition- 
preferably at the 1980 MAGSAT epoch. Effects of non-trivial narrow- 
scale initial magnetic conditions have been considered and tested. The 
truncation of DGRF models is not problematic at non-initial epochs-a 
special feature of the steady flow hypothesis. Such truncation is 
compatible with the goal of matching the evolution of the broad-scale 
radial field (or low-degree scalar potential) . Non-zero damping and a 
velocity field truncation level of 16 were selected. 

Two kinds of weights were investigated: radial field weights and 
general weights. For radial field weights, the weighted, optionally 
constrained, optionally damped iterated linear 1 eas t -squares method 
yields a steady fluid velocity that tracks the weighted evolution of the 
radial magnetic field at Earth's surface indicated by the geomagnetic 
field models fitted. For general weights, the method estimates the 
fluid motion that tracks the weighted evolution of the scalar 
geomagnetic potential at Earth's surface indicated by the geomagnetic 
field models fitted. Both radial field and general weight matrices were 
derived from the covariance matrices for the 1945, 1950, 1955, and 1960 
DGRF models and from the candidate DGRF models for the 1965, 1970, and 
1975 DGRF models. The 1980 DGRF model is treated as comparatively 
perfect-a fair approximation. 

The secular change of the broad-scale radial field relative to 1980 
called for by the DGRF models was analyzed. The results for radial 
field weights are summarized in Table 1. These weights nearly double 
the apparent signal magnitude. The semi -normal i zed total weighted 
variance in the secular change of the broad-scale radial field is 
12,217.8; the tolerable semi -normal ized square-weighted residual is 7. 
Some 99.94271% of the total weighted variance is thus regarded as signal 
to be fit. Normalized root -mean-square residuals of 2.394% or less are 
tolerable for radial field weights. 

The results for general weights are summarized in Table 2. Then 
the semi-normalized total weighted variance in the secular change of the 
low-degree scalar geomagnetic potential relative to 1980 called for by 
the DGRF models radial field is 57,548.6; but the tolerable semi- 
normalized total square-weighted residual is still 7. Some 99.98784% of 
the total generally weighted variance is regarded as signal to be fit. 
Normalized root -mean-square residuals of 1.103% or less are tolerable 
for the generalized weights. 

Both radial field and general weights yield estimated tolerable 
residuals less than the 7% errors expected in the SFM/FFC approximation. 
When using radial field weights, the 7% expected error from this 
approximation will cause a typical misfit of about 2.92a' per DGRF 
epoch. When using general weights, the 7% expected error will cause a 
typical misfit of 6.35a per DGRF epoch. Therefore, the combined 
supposition of a SFM surrounding a FFC in surficially steady flow is 
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expected to be rejectable with appreciable confidence. If this is so, 
the need for a superior earth model will be clear; steps to develop such 
a model have been taken. 

If application of the method developed in this paper to derive a 
single steady flow from 1980 back to 1945 yields errors exceeding 7%, 
then the steady flow hypothesis might reasonably be rejected. If the 
steady flow hypothesis yields errors of 7% or less, then it should be 
retested using a superior earth model. In light of the interim upper 
bound of 7% errors (i.e., at least 99.51% weighted variance reductions) 
on the steady flow hypothesis, the prospects for resolving time- 
dependent flow with an improved earth model by, say, seeking piecewise 
steady flows, may not be too good. The ratio of the time-dependent flow 
signal strength relative to the noise in the definitive secular change 
models would then be at most 2.92 per DGRF epoch for radial field 
weights. But corrections for CMB ellipticity and topography or mantle 
conductivity might reduce the residuals resulting from the supposition 
of steady flow, and, more importantly, this ratio is 6.35 per DGRF epoch 
for general weights. The latter figure seems large enough to hope for 
resolution of time-dependent core flow with a superior earth model. 
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APPENDIX 


* For epochs t n = 1945, 1950, 1955, and 1960, calculation of i\(t ) 

= £2 ij(t n ) was possible for the DGRF models using the derived E- (t ). 
This was not possible for epochs 1965, 1970, and 1975. This might’ve 

due to roundoff error in the published Gauss coefficients, the 
interpretation of the DGRF weighting procedure used to obtain the 
F i j ( 1965-1 97 5 ) , or failure of these E — to account for properly expected 
errors (e.g., fortuitous agreement of candidate models introducing a 
singularity as described in section 3.1). On the chance that random 
roundoff error was the source of this problem, (0.5 nT) 2 was added to 
the diagonals of V matrices for 1965, 1970, and 1975. The resulting V* 

matrices were used to calculate the generalized weight matrices 

«*<t n ) =E*(t n ) _1 = (E(t n ) + R[ (0.5) 2 I]R) _1 

= (R(V(t n ) + 0 . 251] R] -1 = {RvNt^Rr 1 (Al) 

for epochs t ^ = 1965, 1970, and 1975. This procedure also eliminates 

any singularity caused by agreement of the candidate models. The Cl* 
matrices are used to calculate the generalized weighted residual 
variance. The arithmetic mean of these matrices at adjacent epochs is 
inverted to obtain the generalized mean weight matrices Q*(t^) which 
replace W ( t x ) in equations (9) when the generalized weights are used. 

Use of (Al) assumes additional uncorrelated noise in the large DGRF 
models at epochs 1965, 1970, and 1975. The expected amplitude of this 
noise in the radial component 

10 (n+1 ) 2 n 1/2 10 1/2 

( £ £ (0.5 nT) 2 } = { I (n+1 ) 2 ( 0 . 5 nT) 2 } 

n=l 2n+l m=0 n=l 


= 11.24 nT ( rms ) (A2) 

is comparable to that of very long wavelength crustal fields expected by 
VOORHIES (1984) 


1/2 10 1/2 
{E[B rx 2 } } = ( E B rx (n) 2 } 

n=l 


10 1/2 
{ 1 11.91(0.9969)°) 

n=l 


or by LANGEL 
1/2 

{E[B rx 2 } } 


= 10.82 nT 

ESTES & SABAKA (1989), 


10 (n+1) 1/2 

W 


10 (n+1) 


( 2 Rnv) = { I 


n=l 2n+l 


n=l 2n+l 


1/2 

20(0. 9999387 ) n } 


(A3) 


= 10.56 nT. (A4) 

However, the parabolic roundoff spectrum associated with (A2) is much 
more 'blue' than the nearly 'white' crustal spectra (A3) or (A4). 
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